/* -*- mode: C++; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -*- */
// vim:sts=8:sw=8:ts=8:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s

/* Copyright (C) 2010 LinBox
 * Written by <brice.boyer@imag.fr>
 * Adapted from nullspaceMP in Storjohann and Chen's IML.
 *
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 */

/*! @file tests/test-nullspaceZ.C
 * @ingroup tests
 * @brief We test various implementation of dense integer nullsapce.
 * In details :
 *  - ?
 *  - ?
 *  - ?
 */

#define __LINBOX_HAVE_IML

#include "linbox/linbox-config.h"
#include "linbox/field/modular.h"
#include "linbox/integer.h"
#include "linbox/field/PID-integer.h"
#include "linbox/algorithms/matrix-hom.h"
#include "linbox/algorithms/blas-domain.h"
#include "linbox/matrix/random-matrix.h"
#include "linbox/randiter/generic.h"
#include "linbox/algorithms/NullspaceZ.h"
#include "linbox/util/timer.h"

#include "iml.h"

using namespace LinBox;

#if 0
template<class BlackBox>
bool isZeroProductRandom(const BlackBox & A, const BlackBox & B)
{
	size_t m = A.rowdim();
	size_t k = A.coldim();
	size_t n = B.coldim();

	linbox_check( k == B.rowdim() );
}
#endif

/*! @todo si corps fini > on teste pour différent vecteurs
 *  si entiers > on teste pour suffisament de nombre premiers
 *  si blackbox, on est content
 *  si non, on fait comme si
 */

// n^{2} time random test
template<class _Field>
bool isProductZeroRandom(const DenseMatrix<_Field> & A, const DenseMatrix<_Field> & B)
{
	size_t m = A.rowdim();
	size_t k = A.coldim();
	size_t n = B.coldim();

	linbox_check( k == B.rowdim() );

	typedef Modular<double> Field ;
	size_t bound = 3 ; //! @bug  dépend de la taille max
	std::vector<double> x(n);
	std::vector<double> y(k);
	std::vector<double> z(m);
	for (size_t i = 0 ; i < bound ; ++i)
	{
		Integer p = random() ;
		Field F(p);
		BlasMatrix<typename Field::Element> Ap (m,k);
		BlasMatrix<typename Field::Element> Bp (k,n);
		// DenseMatrix<Field> Ap (F,m,k);
		// DenseMatrix<Field> Bp (F,k,n);

		MatrixHom::map(Ap,A,F);
		MatrixHom::map(Bp,B,F);
		VectorDomain<Field> VD(F);
		VD.random(x);
		BlasMatrixDomain<Field> MD(F);
		// MatrixDomain<Field> BMD(F);
		// BMD.vectorMul(y,Bp,x) ;
		// BMD.vectorMul(z,Ap,y) ;
		MD.mul(y,Bp,x);
		MD.mul(z,Ap,y);
		/*   Ap.apply(z,y); */ ///@todo  whyyyyy ?
		if (! VD.isZero(z) )
			return false ;
	}
	return true;
}

template<class Element>
bool isProductZeroRandom(const BlasMatrix<Element> & A, const BlasMatrix<Element> & B)
{
	size_t m = A.rowdim();
	size_t k = A.coldim();
	size_t n = B.coldim();

	linbox_check( k == B.rowdim() );

	typedef Modular<double> Field ;
	size_t bound = 3 ; //! @bug  dépend de la taille max
	std::vector<double> x(n);
	std::vector<double> y(k);
	std::vector<double> z(m);
	for (size_t i = 0 ; i < bound ; ++i)
	{
		Integer p = Integer::random_between(15,18) ;
		Field F(p);
		BlasBlackbox<Field> Ap (F,m,k);
		BlasBlackbox<Field> Bp (F,k,n);

		MatrixHom::map(Ap,A,F);
		MatrixHom::map(Bp,B,F);
		VectorDomain<Field> VD(F);
		VD.random(x);
		Bp.apply(y,x); // Compose ?
		Ap.apply(z,y);
		if (! VD.isZero(z) ){
			std::cout << "\033[1;31miznogood\033[0m" << std::endl;
			return false ;
		}
	}
	std::cout << "\033[1;32mis cool\033[0m" << std::endl;
	return true;
}

template<class Element>
bool isProductZeroRandom(const BlasMatrix<Element> & A, const BlasMatrix<Element> & B,
			 double p, int d)
{
	size_t m = A.rowdim();
	size_t k = A.coldim();
	size_t n = B.coldim();

	linbox_check( k == B.rowdim() );

	typedef Modular<Integer> Field ;
	size_t bound = 3 ; //! @bug  dépend de la taille max
	std::vector<Integer> x(n);
	std::vector<Integer> y(k);
	std::vector<Integer> z(m);
	Integer pd = pow((integer)p,(long)d) ;
	Field F(pd);
	for (size_t i = 0 ; i < bound ; ++i)
	{
		BlasBlackbox<Field> Ap (F,m,k);
		BlasBlackbox<Field> Bp (F,k,n);

		MatrixHom::map(Ap,A,F);
		MatrixHom::map(Bp,B,F);
		VectorDomain<Field> VD(F);
		VD.random(x);
		Bp.apply(y,x); // Compose ?
		Ap.apply(z,y);
		if (! VD.isZero(z) ){
			std::cout << "\033[1;31miznogood\033[0m" << std::endl;
			return false ;
		}
	}
	std::cout << "\033[1;32mis cool\033[0m" << std::endl;
	return true;
}


//! todo : BlasBlackbox(PID_integer,
template<class BlackBox>
bool isProductZeroRandom(const BlackBox & A, const BlackBox & B)
{
	size_t m = A.rowdim();
	size_t k = A.coldim();
	size_t n = B.coldim();

	linbox_check( k == B.rowdim() );

	typedef Modular<double> Field ;
	size_t bound = 3 ; //! @bug  dépend de la taille max
	std::vector<double> x(n);
	std::vector<double> y(k);
	std::vector<double> z(m);
	for (size_t i = 0 ; i < bound ; ++i)
	{
		Integer p = random() ;
		Field F(p);
		BlasBlackbox<Field> Ap (m,k);
		BlasBlackbox<Field> Bp (k,n);

		MatrixHom::map(Ap,A,F);
		MatrixHom::map(Bp,B,F);
		VectorDomain<Field> VD(F);
		VD.random(x);
		Bp.apply(y,x);
		Ap.apply(z,y);
		if (! VD.isZero(z) )
			return false ;
	}
	return true;
}


/** @brief tests integer nullspace routines.
 *
 * @param m rows
 * @param n cols
 * @param l size of random integers (\f$0,2^l\f$)
 */
int
test_NullspaceIML( const size_t m, const size_t n, const size_t l, size_t r)
{
	Timer tim ;
	//BlasMatrix<Integer> A(m,n);
	PID_integer IntRing ;
	// typedef DenseMatrixBase<PID_integer> IntMat ;
	// IntMat A(IntRing,m,n);
	typedef BlasMatrix<Integer> IntMat ;
	IntMat A(m,n);
	// int r = random() %(std::min(m,n)) ;
	RandomIntegerIter<> R(l);
	// typedef PID_integer::RandIter R ;
	typedef RandomIntegerIter<> IntRandIter ;

	typedef RandomDenseMatrix<IntRandIter, PID_integer > IntRand_t;

	IntRand_t Randomize(IntRing,R);
	// Randomize.random<IntMat>(A);
	// A.write(std::cout) << std::endl;
	//

	// size_t r = Integer::random_lessthan((Integer)std::min(m,n)) ;
	// size_t r = std::min(m,n)-1 ;

	// std::cout << "probable target rank : " << r << std::endl;

	// Randomize.randomRank<IntMat>(A,r,RankBuilder::Rank_update_());
	Randomize.randomRank<IntMat>(A,r,RankBuilder::LU_());
	// A.write(std::cout << "Input matrix\n A:=" << std::endl,true) << ';' << std::endl;
	// return EXIT_FAILURE ;
	size_t lA ;
	MaxElementSize(lA,A);
	std::cout << "A max log-size: " << lA << std::endl;
	// if (lA != l) std::cout << " *** WARNING *** expected size : " << l << " but got " << lA << std::endl;

	/*  IML copy of A */
#ifdef __LINBOX_HAVE_IML
	BlasMatrix<Integer> AA((const BlasMatrix<Integer>)A);
	mpz_t * AA_ptr = reinterpret_cast<mpz_t*>(AA.getPointer());
	// for (size_t i = 0 ; i < m ; ++i) {
	// for (size_t j = 0 ; j < n ; ++j) {
	// Integer Aij ;
	// mpz_set(Aij.get_mpz(),AA_ptr[i*n+j]);
	// if (Aij != AA.getEntry(i,j))
	// std::cout << 'x' ;
	// else
	// std::cout << 'o' ;
	// }
	// std::cout << std::endl;
	// }
#endif

#if 0
	mpz_t * AA_ptr = new mpz_t[m*n] ;
	for (size_t i = 0 ; i < m ; ++i)
		for (size_t j = 0 ; j < n ; ++j)
			mpz_init_set(AA_ptr[i*n+j],
				     (mpz_ptr)(LinBox::SpyInteger::get_rep(A.getEntry(i,j))));
#endif
	/*  ------------ */


	BlasMatrix<Integer> Z ;
	// DenseMatrix<PID_integer> Z(IntegerRing) ;

	// std::vector<Integer> d;
	tim.clear() ; tim.start();
	size_t t = IntegerNS(Z/*,d*/,A,LinBoxTag::Right);
	tim.stop(); std::cout << "\033[0;33m" << tim << "\033[0m" << std::endl;

	SumElementSize(lA,Z);
	std::cout << " Z of log-sum-size : " << ((lA)?(/*std::log*/(lA)):(0)) << std::endl;
	std::cout << " *** ker dim " << t  << " ***" << std::endl;
	// right
	bool null ;
	if (!Z.coldim()) {
		if (t != std::min(m,n)) {
			std::cout << "\033[1;31miznogood\033[0m" << std::endl;
			null = false ;
		}
		else
			null = true;
	}
	else
		null = isProductZeroRandom(A,Z);

	if (!null)
		return EXIT_FAILURE ;

#ifdef __LINBOX_HAVE_IML
	std::cout << ">>IML<<" << std::endl;
	mpz_t ** ZZ_ptr = new mpz_t* ;

	Timer er ;
	er.clear() ; er.start();
	size_t s = nullspaceMP(m,n,AA_ptr,ZZ_ptr);
	er.stop(); std::cout << "\033[0;33m" << er << "\033[0m" << std::endl;

	BlasMatrix<Integer> ZZ(n,s);
	for (size_t i = 0 ; i < n ; ++i)
		for (size_t j = 0 ; j < s ; ++j) {
			Integer Zij ;
			mpz_set(Zij.get_mpz(),(*ZZ_ptr)[i*s+j]);
			ZZ.setEntry(i,j,Zij);
		}
	size_t lB ;
	SumElementSize(lB,ZZ);
	std::cout << "ZZ of log-sum-size : " << ((lB)?(/*std::log*/(lB)):(0)) << std::endl;
	std::cout << "gain : \033[1;34m" << double(lB-lA)/lB*100. << "\%\033[0m over IML (in bits)" << std::endl;
	std::cout << " *** ker dim " << s  << " ***" << std::endl;

	if (!ZZ.coldim()) {
		if (s != std::min(m,n)) {
			std::cout << "\033[1;31miznogood\033[0m" << std::endl;
			null = false ;
		}
		else
			null = true;
	}
	else
		null = isProductZeroRandom(A,ZZ);

	if (!null)
		return EXIT_FAILURE ;

	// for (size_t i = 0 ; i < m*n ; ++i)
	// mpz_clear(AA_ptr[i]);
	// delete[] AA_ptr ;
	for (size_t i = 0 ; i < s*n ; ++i)
		mpz_clear((*ZZ_ptr)[i]);
	delete ZZ_ptr ;

	if (tim.realtime() > er.realtime())
		std::cout << "\033[0;31m *** IML is " << tim.realtime()/er.realtime() << " faster than LinBox\033[0m" << std::endl;
	else
		std::cout << "\033[0;32m *** LinBox is " << er.realtime()/tim.realtime() << " faster than IML\033[0m" << std::endl;
#endif
	return EXIT_SUCCESS ;
}

int
test_NullspaceCRA( const size_t m, const size_t n, const size_t l, size_t r)
{
	std::cout << "cra" << std::endl;
	Timer tim ;
	//BlasMatrix<Integer> A(m,n);
	PID_integer IntRing ;
	// typedef DenseMatrixBase<PID_integer> IntMat ;
	// IntMat A(IntRing,m,n);
	typedef BlasMatrix<Integer> IntMat ;
	IntMat A(m,n);
	// int r = random() %(std::min(m,n)) ;
	RandomIntegerIter<> R(l);
	// typedef PID_integer::RandIter R ;
	typedef RandomIntegerIter<> IntRandIter ;

	typedef RandomDenseMatrix<IntRandIter, PID_integer > IntRand_t;

	IntRand_t Randomize(IntRing,R);
	// Randomize.random<IntMat>(A);
	// A.write(std::cout) << std::endl;
	//

	// size_t r = Integer::random_lessthan((Integer)std::min(m,n)) ;
	// size_t r = std::min(m,n)-1 ;

	// std::cout << "probable target rank : " << r << std::endl;

	// Randomize.randomRank<IntMat>(A,r,RankBuilder::Rank_update_());
	Randomize.randomRank<IntMat>(A,r,RankBuilder::LU_());
	// A.write(std::cout << "Input matrix\n A:=" << std::endl,true) << ';' << std::endl;
	// return EXIT_FAILURE ;
	size_t lA ;
	MaxElementSize(lA,A);
	std::cout << "A max log-size: " << lA << std::endl;
	// if (lA != l) std::cout << " *** WARNING *** expected size : " << l << " but got " << lA << std::endl;

	/*  IML copy of A */
#ifdef __LINBOX_HAVE_IML
	BlasMatrix<Integer> AA((const BlasMatrix<Integer>)A);
	mpz_t * AA_ptr = reinterpret_cast<mpz_t*>(AA.getPointer());
	// for (size_t i = 0 ; i < m ; ++i) {
	// for (size_t j = 0 ; j < n ; ++j) {
	// Integer Aij ;
	// mpz_set(Aij.get_mpz(),AA_ptr[i*n+j]);
	// if (Aij != AA.getEntry(i,j))
	// std::cout << 'x' ;
	// else
	// std::cout << 'o' ;
	// }
	// std::cout << std::endl;
	// }
#endif

#if 0
	mpz_t * AA_ptr = new mpz_t[m*n] ;
	for (size_t i = 0 ; i < m ; ++i)
		for (size_t j = 0 ; j < n ; ++j)
			mpz_init_set(AA_ptr[i*n+j],
				     (mpz_ptr)(LinBox::SpyInteger::get_rep(A.getEntry(i,j))));
#endif
	/*  ------------ */


	BlasMatrix<Integer> Z ;
	// DenseMatrix<PID_integer> Z(IntegerRing) ;

	// std::vector<Integer> d;
	tim.clear() ; tim.start();
	size_t t = RightCradicNullspace<IntMat >(Z/*,d*/,A,Dixon_);
	tim.stop(); std::cout << "\033[0;33m" << tim << "\033[0m" << std::endl;
	Z.write(std::cout << "Z :=",true) << ';' << std::endl;
	A.write(std::cout << "A :=",true) << ';' << std::endl;

	SumElementSize(lA,Z);
	std::cout << " Z of log-sum-size : " << ((lA)?(/*std::log*/(lA)):(0)) << std::endl;
	std::cout << " *** ker dim " << t  << " ***" << std::endl;
	// right
	bool null ;
	if (!Z.coldim() || !Z.rowdim()) {
		if (t != std::min(m,n)) {
			std::cout << "\033[1;31miznogood\033[0m" << std::endl;
			null = false ;
		}
		else
			null = true;
	}
	else
		null = isProductZeroRandom(A,Z);

	if (!null)
		return EXIT_FAILURE ;

#ifdef __LINBOX_HAVE_IML
	std::cout << ">>IML<<" << std::endl;
	mpz_t ** ZZ_ptr = new mpz_t* ;

	Timer er ;
	er.clear() ; er.start();
	size_t s = nullspaceMP(m,n,AA_ptr,ZZ_ptr);
	er.stop(); std::cout << "\033[0;33m" << er << "\033[0m" << std::endl;

	BlasMatrix<Integer> ZZ(n,s);
	for (size_t i = 0 ; i < n ; ++i)
		for (size_t j = 0 ; j < s ; ++j) {
			Integer Zij ;
			mpz_set(Zij.get_mpz(),(*ZZ_ptr)[i*s+j]);
			ZZ.setEntry(i,j,Zij);
		}
	size_t lB ;
	SumElementSize(lB,ZZ);
	std::cout << "ZZ of log-sum-size : " << ((lB)?(/*std::log*/(lB)):(0)) << std::endl;
	std::cout << "gain : \033[1;34m" << double(lB-lA)/lB*100. << "\%\033[0m over IML (in bits)" << std::endl;
	std::cout << " *** ker dim " << s  << " ***" << std::endl;

	if (!ZZ.coldim()) {
		if (s != std::min(m,n)) {
			std::cout << "\033[1;31miznogood\033[0m" << std::endl;
			null = false ;
		}
		else
			null = true;
	}
	else
		null = isProductZeroRandom(A,ZZ);

	if (!null)
		return EXIT_FAILURE ;

	// for (size_t i = 0 ; i < m*n ; ++i)
	// mpz_clear(AA_ptr[i]);
	// delete[] AA_ptr ;
	for (size_t i = 0 ; i < s*n ; ++i)
		mpz_clear((*ZZ_ptr)[i]);
	delete ZZ_ptr ;

	if (tim.realtime() > er.realtime())
		std::cout << "\033[0;31m *** IML is " << tim.realtime()/er.realtime() << " faster than LinBox\033[0m" << std::endl;
	else
		std::cout << "\033[0;32m *** LinBox is " << er.realtime()/tim.realtime() << " faster than IML\033[0m" << std::endl;
#endif
	return EXIT_SUCCESS ;
}


int
test_Nullspace_Padic( const size_t m, const size_t n, const size_t l, double p, double d, size_t r,
		      NullSpaceMethod meth)
{
	Timer tim ;

	//BlasMatrix<Integer> A(m,n);
	PID_integer IntRing ;
	// typedef DenseMatrixBase<PID_integer> IntMat ;
	// IntMat A(IntRing,m,n);
	typedef BlasMatrix<Integer> IntMat ;
	IntMat A(m,n);
	// int r = random() %(std::min(m,n)) ;
	RandomIntegerIter<> R(l);
	// typedef PID_integer::RandIter R ;
	typedef RandomIntegerIter<> IntRandIter ;

	typedef RandomDenseMatrix<IntRandIter, PID_integer > IntRand_t;

	IntRand_t Randomize(IntRing,R);
	// Randomize.random<IntMat>(A);
	// A.write(std::cout) << std::endl;
	//

	// size_t r = Integer::random_lessthan((Integer)std::min(m,n)) ;

	// std::cout << "probable target rank : " << r << std::endl;

	// Randomize.randomRank<IntMat>(A,r,RankBuilder::Rank_update_());
	Randomize.randomRank<IntMat>(A,r,RankBuilder::LU_());
	// A.write(std::cout << "#Input matrix\n A:=",true) << ';' << std::endl;
	// return EXIT_FAILURE ;

	size_t lA ;
	MaxElementSize(lA,A);
	std::cout << "A : " << m << 'x' << n << " (" << lA << ") [" << p << '^' << d << ']' << std::endl;
	// if (lA != l) std::cout << " *** WARNING *** expected size : " << l << " but got " << lA << std::endl;


	BlasMatrix<Integer> Z ;
	// DenseMatrix<PID_integer> Z(IntegerRing) ;

	// std::vector<Integer> d;

	tim.clear() ; tim.start();
	size_t t = RightPadicNullspace (Z/*,d*/,A,p,d,meth);
	tim.stop(); std::cout << "\033[0;33m" << tim << "\033[0m" << std::endl;

	std::cout << " *** ker dim " << t  << " ***" << std::endl;
	// right
	bool null ;
	if (!Z.coldim()){
		if (t != std::min(m,n)) {
			std::cout << "\033[1;31miznogood\033[0m" << std::endl;
			null = false ;
		}
		else
			null = true;
	}
	else
		null = isProductZeroRandom(A,Z,p,d);

	if (!null)
		return EXIT_FAILURE ;

	return EXIT_SUCCESS ;
}
#if 0
int
test_Nullspace_Noix( const size_t m, const size_t n, const size_t l, double p, double d, size_t r)
{
	Timer tim ;

	//BlasMatrix<Integer> A(m,n);
	PID_integer IntRing ;
	// typedef DenseMatrixBase<PID_integer> IntMat ;
	// IntMat A(IntRing,m,n);
	typedef BlasMatrix<Integer> IntMat ;
	IntMat A(m,n);
	// int r = random() %(std::min(m,n)) ;
	RandomIntegerIter<> R(l);
	// typedef PID_integer::RandIter R ;
	typedef RandomIntegerIter<> IntRandIter ;

	typedef RandomDenseMatrix<IntRandIter, PID_integer > IntRand_t;

	IntRand_t Randomize(IntRing,R);
	// Randomize.random<IntMat>(A);
	// A.write(std::cout) << std::endl;
	//

	// size_t r = Integer::random_lessthan((Integer)std::min(m,n)) ;
	// size_t r = std::min(m,n)-1 ; //Integer::random_lessthan((Integer)std::min(m,n)) ;

	// std::cout << "probable target rank : " << r << std::endl;

	// Randomize.randomRank<IntMat>(A,r,RankBuilder::Rank_update_());
	Randomize.randomRank<IntMat>(A,r,RankBuilder::LU_());
	size_t lA ;
	MaxElementSize(lA,A);
	std::cout << "A max log-size: " << lA << std::endl;

	// if (lA != l) std::cout << " *** WARNING *** expected size : " << l << " but got " << lA << std::endl;
	// A.write(std::cout << "Input matrix\n A:=",true) << ';' << std::endl;
	// return EXIT_FAILURE ;



	BlasMatrix<Integer> Z ;
	// DenseMatrix<PID_integer> Z(IntegerRing) ;

	// std::vector<Integer> d;

	tim.clear();tim.start();
	size_t t = RightNoix(Z/*,d*/,A,p,d);
	tim.stop(); std::cout << "\033[0;33m" << tim << "\033[0m" << std::endl;


	std::cout << " *** ker dim " << t  << " ***" << std::endl;
	// right
	bool null ;
	if (!Z.coldim()){
		if (t != std::min(m,n)) {
			std::cout << "\033[1;31miznogood\033[0m" << std::endl;
			null = false ;
		}
		else
			null = true;
	}
	else
		null = isProductZeroRandom(A,Z,p,d);

	if (!null)
		return EXIT_FAILURE ;

	return EXIT_SUCCESS ;
}

int
test_Nullspace_MoenckCarter( const size_t m, const size_t n, const size_t l, double p, double d, size_t r)
{
	Timer tim ;

	//BlasMatrix<Integer> A(m,n);
	PID_integer IntRing ;
	// typedef DenseMatrixBase<PID_integer> IntMat ;
	// IntMat A(IntRing,m,n);
	typedef BlasMatrix<Integer> IntMat ;
	IntMat A(m,n);
	// int r = random() %(std::min(m,n)) ;
	RandomIntegerIter<> R(l);
	// typedef PID_integer::RandIter R ;
	typedef RandomIntegerIter<> IntRandIter ;

	typedef RandomDenseMatrix<IntRandIter, PID_integer > IntRand_t;

	IntRand_t Randomize(IntRing,R);
	// Randomize.random<IntMat>(A);
	// A.write(std::cout) << std::endl;
	//

	// size_t r = Integer::random_lessthan((Integer)std::min(m,n)) ;
	// size_t r = std::min(m,n)-1 ; //Integer::random_lessthan((Integer)std::min(m,n)) ;

	// std::cout << "probable target rank : " << r << std::endl;

	// Randomize.randomRank<IntMat>(A,r,RankBuilder::Rank_update_());
	Randomize.randomRank<IntMat>(A,r,RankBuilder::LU_());
	size_t lA ;
	MaxElementSize(lA,A);
	std::cout << "A max log-size: " << lA << std::endl;

	// if (lA != l) std::cout << " *** WARNING *** expected size : " << l << " but got " << lA << std::endl;
	// A.write(std::cout << "Input matrix\n A:=",true) << ';' << std::endl;
	// return EXIT_FAILURE ;



	BlasMatrix<Integer> Z ;
	// DenseMatrix<PID_integer> Z(IntegerRing) ;

	// std::vector<Integer> d;

	tim.clear();tim.start();
	size_t t = RightMoenckCarter(Z/*,d*/,A,p,d);
	tim.stop(); std::cout << "\033[0;33m" << tim << "\033[0m" << std::endl;

	std::cout << " *** ker dim " << t  << " ***" << std::endl;
	// right
	bool null ;
	if (!Z.coldim()){
		if (t != std::min(m,n)) {
			std::cout << "\033[1;31miznogood\033[0m" << std::endl;
			null = false ;
		}
		else
			null = true;
	}
	else
		null = isProductZeroRandom(A,Z,p,d);

	if (!null)
		return EXIT_FAILURE ;

	return EXIT_SUCCESS ;
}
#endif
int lauch_test(NullSpaceMethod Meth)
{
	float prop = 0.6 ;
	int res = 0 ;
	int n = 100 ;
	int l = 100 ;
	int r = prop*n ;
	// int p, int d;

#if 0 /*  TESTS */
	n = 10, r=prop*n , l = 10 ;
	res |= test_Nullspace_Padic(n,n,l,11,3,r,Meth);
	res |= test_Nullspace_Padic(n,n,3*l,11,10,r,Meth);
	res |= test_Nullspace_Padic(n,n,10*l,11,10,r,Meth);
	// n = 500 ; r = prop*n ; l = 5000 ;
	// res |= test_Nullspace_Padic(n,n,l,65521,5,r,Meth);
	// res |= test_Nullspace_Padic(n,n,l,65521,10,r,Meth);
	// res |= test_Nullspace_Padic(n,n,l,65521,30,r,Meth);
#endif

#if 1 /*  BENCH */
	/*  */
	n = 500 ; r = prop*n ; l = 500 ;
	res |= test_Nullspace_Padic(n,n,l, 65521,5 ,r,Meth);
	res |= test_Nullspace_Padic(n,n,l, 65521,10,r,Meth);
	res |= test_Nullspace_Padic(n,n,l, 65521,20,r,Meth);
	/*   */
	n = 800 ; r = prop*n ; l = 1000 ;
	res |= test_Nullspace_Padic(n,n,l, 65521,5 ,r,Meth);
	res |= test_Nullspace_Padic(n,n,l, 65521,10,r,Meth);
	res |= test_Nullspace_Padic(n,n,l ,65521,20,r,Meth);
	/*   */
	n = 1500 ; r = prop*n ;
	res |= test_Nullspace_Padic(n,n,l, 65521,5 ,r,Meth);
	res |= test_Nullspace_Padic(n,n,l, 65521,10,r,Meth);
	res |= test_Nullspace_Padic(n,n,l, 65521,20,r,Meth);
	/*   */
	n = 2000 ; r = prop*n ;
	res |= test_Nullspace_Padic(n,n,l, 65521,5 ,r,Meth);
	res |= test_Nullspace_Padic(n,n,l, 65521,10,r,Meth);
	res |= test_Nullspace_Padic(n,n,l, 65521,20,r,Meth);
	// n = 5000 ; r = prop*n ; l = 500 ;
	// res |= test_Nullspace_Padic(n,n,l,65521,5,r,Meth);
	// res |= test_Nullspace_Padic(n,n,l,65521,30,r,Meth);
#endif
	return res ;

}
int main()
{
	srand(time(NULL));
	srand48(time(NULL));
	size_t _size = 10 ;
	size_t n = _size ;
	size_t r = 0.6*n  ;
	linbox_check(r<=n);
	// size_t r = Integer::random_lessthan((Integer)std::min(m,n)) ;
	// size_t r = std::min(m,n)-1 ; //Integer::random_lessthan((Integer)std::min(m,n)) ;

	int titi=0,toto=0,tata=0,tutu=0 ;
	std::cout << "\033[1;35m>>>>>> Standard <<<<<<\033[0m" << std::endl;
	// toto = test_NullspaceIML(n,n,10,r);
	// toto = test_NullspaceIML(n,n,40,r);
	// toto = test_NullspaceIML(n,n,100,r);
	// toto = test_Nullspace(n,n,200,r);

	// toto = test_NullspaceCRA(n,n,10,r);
	// toto = test_NullspaceCRA(n,n,40,r);
	// toto = test_NullspaceCRA(n,n,100,r);
#if 1
	std::cout << "\033[1;35m>>>>>> Dixon <<<<<<\033[0m" << std::endl;
	titi = lauch_test(Dixon_);
	std::cout << "\033[1;35m>>>>>> MoenckCarter <<<<<<\033[0m" << std::endl;
	tutu = lauch_test(MoenckCarter_);
	std::cout << "\033[1;35m>>>>>> Noix <<<<<<\033[0m" << std::endl;
	tata = lauch_test(Noix_);
#endif
	return toto||titi||tata||tutu ;
}
